

options(error=stop)
rm(list=ls())

source("1.0 MyFunc.R");  source("1.1.1 TrueValue.R")
source("1.0 Delta_est.R");  source("1.1 DataGen.R")
source("1.2 Reg.R"); source("1.3 IPW.R"); source("1.4 G-est.R")
source("1.5 TR.R")
library(brm)

REP = 1

paras = list(
    rep.cand = 1:1000,
    pi_true = c(1,2),       # 1 = correct, 2 = incorrect
    deltad_true = c(1,2),   # 1 = correct, 2 = incorrect
    delta_true = c(1,2),    # 1 = correct, 2 = incorrect
    op_true = c(1,2),       # 1 = correct, 2 = incorrect
    n.cand = c(500)
)
para.values = ParaValues(paras)

### Comparison method
#  1 = Regression
#  2 = IPW
#  3 = Bounded IPW
#  4 = G-estimation
#  5 = Triply Robust
#  6 = Bounded triply robust
method   <<- c("reg","ipw","b-ipw","g","tr","b-tr")

MESSAGE <-- FALSE     # compress message (This is a global variable!)

# temp = list(output=NULL, params.x = NULL)

ptm0 = ptm1 = proc.time()
result = apply(para.values, 1:length(paras), function(para){
    cat(Message(paras, para, ptm1))
     
    ### 0: Data generation
    data = DataGen(para["n"],SEED=REP * 100 + para["rep"])
    MyAttach(data)
    
    ### 1: Estimation
    Delta.est.sol = DeltaEst(x.delta = xs[[para["delta_true"]]],
            x.delta.d = xs[[para["deltad_true"]]], x.op.y = xs[[para["op_true"]]],
            x.op.d = xs[[para["op_true"]]], x.pi = xs[[para["pi_true"]]], z, d, y)
    return(Delta.est.sol$Delta.est)                                 

})
proc.time() - ptm1
proc.time() - ptm0
dimnames(result) = c(list(method),NameList(paras))
names(dimnames(result)) = c("estimator", NameFirst(names(paras)))
d = length(dim(result))
result = aperm(result, c(2:d,1))    # put method at the last dimension


ptm = proc.time()


save(result, Delta.true, file=paste("./RDataFiles/results",REP,
                                    ".RData", sep = ""))
